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ABSTRACT 



The solution of the nonlinear two-dimensional reactor dynamics 

equation subjected to prompt feedback conditions using the finite 

element technique leads to the matrix formulation A. , \b , = B # .iK + 

M ij j ij J 

(i,j,k = 1, N). This system has been solved directly 

ijk j k 

in a previous work; but because the nonlinearity C ijk^j^k was 
premultiplied by [A] \ large computer storage was required for the 
small problem considered. The task of this thesis is the development 
of computational techniques which allow the problem to be solved for 
large systems. Specifically, these techniques are: (1) the treatment 

of the nonlinearity on the element level, (2) the compacting of the 
sparce matrices to include only non-zero terms, and (3) the construction 
of a new computer code based on the Crank-Nicolson formulation for the 
solution of differential equations. 

To support the theory presented, test problems were solved by the 
original method, the linearized technique, and the Crank-Nicolson 
treatment. The results were analyzed and compared graphically. All 
three of the innovations developed in this thesis appear to be useful 
tools for solving nonlinear time dependent differential equations. 
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I. INTRODUCTION 



When temperature-dependent feedback is considered in the nuclear 
reactor dynamics problem, a nonlinear field equation in space and time 
results [Ref. 1]. For the non-homo geneous , or multi-region reactor, 
however, the space-dependent dynamics behavior following a nonlinear 
initial disturbance is no longer reachable in analytical form. 
Fortunately, by modeling the reactor as a system of finite regions of 
interest where the neutronic properties are known, the method of finite 
elements can be applied to yield the solutions. The fundamental concepts 
relating reactor behavior to the finite element formulation have been 
presented in 1974 [Ref. 2], but a more recent work by Nguyen and 
Salinas [Ref. 3] gives a thorough discussion of the complete problem. 

That work forms the basis and starting point for this thesis, the 
objective of which is to develop improved computational methods for 
dealing with that finite element formulation. 

In Reference 3 the dynamic flux equation under prompt feedback 
conditions was given as 



. dip (r,z,t) p „ 

77 — 2- « D V ij< + Y Z 'p - a K Z ip 

V 3 1 m m m am m mmamm 

m 



a) 



for each non-homogeneous zone "m" contained in the reactor body, 
where the usual symbols are employed: 



X = v /Z - 1 
m fm am 



K = e /p C (°C/unit flux-sec) 
m tm m pm 



Y = (A/V ) <h /p C ) sec 
m m m m m pro 



-i 
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p = heat capacity 

a = reactivity temperature coefficient °C *** 

h = convection heat transfer coefficient 

A/V = heat transfer area to volume of energy 
generation ratio 

= neutron fission cross-section 

= neutron absorption cross-section 

v = neutrons emitted per fission 

e = energy produced per fission 

The subscript "m" will be omitted from further discussion for simplicity. 
The transformation of this general equation into a problem in 

finite dimensional vector space begins by constructing the following N 

term approximation'* - 

N 

I|)(x,t) = <Kx,t) = 2L ^.(t)G (x) (2) 

3=1 J J 

where N is the number of degrees of freedom, or nodes, in the space, 
and the G ^ are the basis functions of the approximate solution space . 
The Galerkin method seeks to make the residual 

R(x,t) = L ^ - f, 

where L\f> = f is the field equation, othogonal to each of the basis 
functions so that 



This discussion leading to the establishment of the matrix 
equations is essentially an abridgement of the development given in 
reference 3. 
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/ G (x)R(x,t) 
J S 



dx = 0 



i = 1, 2, 



N . 



(3) 



For the nuclear reactor dynamics problem of Eq. (1), the residual 
becomes 



R(2L, t) 



li 2 

at 



- V DV 2 <j> - VXZ i + wZ i 2 
a a 



(4) 



with VaKE = w . Putting this expression into Eq. (3) and integrating 
a 

by parts (a distinct advantage of Galerkin) , the following coefficients 
emerge 



A IJ ■ ff G I G K rdrdz 



(5a) 



IJ 



-ff- 

S L 



3G. 

3r" 



3r 



+ 



3G. 

3z 



3G^ 

az 



rdrdz 



(5b) 



IJK 



■ff 



g i g j g k 



rdrdz 



(5c) 



I.J.K == 1, 



N 



in cylindrical coordinates where dS = rdrdz. This allows Eq. (1) to be 
written as follows for the case of uniform neutronic properties within each 
reactor region m. 



N N N 

Z a iA - - V I D I Z B u*j + Yi 1 .! Z A u*j 

J=1 J=1 J=1 



N N 

z z 

J=1 K=1 



“i Z Z c uk*A I = 1 > 2 ' ••• fI 



( 6 ) 
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This becomes, upon combining constants, in Einstein summation convention, 
where there is no sum on an underlined repeated index, 



a u *j ■ ab u ♦j - “ 




(7a) 






I, J,K, = 1 N 



With boundary conditions, the system is now well posed and ready for 

solution. The equation-solver must be chosen with care, however, 

because this system of ordinary differential equations is both stiff 

and non-linear. "Stiff" is used here and throughout this thesis to 

denote a system which gives a large response to a small stimulus; in 

this case, the nature of the reactor dynamics problem predicts a large 

2 

change in flux over a small change in time. Current experience 
suggests Gear’s predictor-corrector method [Ref. 4] as written in the 
computer programming code DVOGER [Ref. 5]. It requires only 1) the 
locus (value) of the points at a given time, 2) a routine to evaluate 
the instantaneous derivatives at any time and 3) a routine to evaluate 
the Jacobian (J = 3f^(iJj, t) /3t/; ) . Based on this information, the program 
assumes a time step, which may be quite small (the minimum size being a 
function of the particular computer) , and tries to fit a trial (predicted) 
solution for that time interval. The corrected (integrated) solution is 
then obtained by iteration until convergence is attained. If the 
predicted and corrected solutions fail to match within a specified error 
criterion, the time step is decreased and the process repeated. On the 
other hand, if "excessive" accuracy is attained, the next attempted 
time step will be increased. Gear’s method in DVOGER is well suited 
to stiff non-linear systems. Also contained in DVOGER is the Adam’s 

2 

This definition is somewhat broader in scope than that used by 
many other authors. 
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method which parallels that of Gear except that the Jacobian is not 
computed. Lacking the additional information about the rate of change, 
the Adam’s method must proceed much more cautiously with stiff systems 
and hence marches forward with exceedingly small time steps. Present 
experience confirms the caveat advice given by DVOGER that the Adam’s 
method is not intended for ’’stiff” cases. 

The difficulty that arises from acquiring solutions to Eq. (7), 
then, should not, and, in light of the present experience, does not 
result from DVOGER or any other acceptable equation solver, but is, rather, 
a function of matrix size, vis a vis the number of nodal points in the 
finite element approximation. The finite element modeling of large 
reactors or the close scrutiny of small ones is severely restricted, 
therefore, unless the number of nodal points can be raised significantly. 

As an example, the test case considered in this thesis consisted of only 
38 nodes but, \vdien processed directly from Eq. (7), required over 300K 
bytes of storage in the IBM 360/67 computer. This requirement came 
largely from the [A] ^*[C ] matrix alone needing a size of 4(38x38x38) = 

219K bytes in single precision. Manipulations with this cubic were 
indeed quite burdensome. Since machine processing time is, other things 
being equal, dependent on size, the combined time and space requirements 
would prohibit the consideration of even moderate sized problems except 
on the largest computers. Additionally, the inversion of [A] in Eq. 

(7) is indicated in order to provide an explicit value of ip for the 
DVOGER routine. It is this difficulty of size, and with it processing 
time, to which the remainder of this thesis is devoted. Experience with 
the DVOGER routine is presented, and the author's equation-solver based on 
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a Crank-Nicolson formulation is discussed as the analysis of the size 
problem is developed. Table I shows clearly the influence of matrix 
size on computer storage requirements. 
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II. REDUCTION OF THE NONLINEARITY 



A. GLOBAL TREATMENT OF THE NONLINEARITY 

The evaluation of Eq. (7) generates the following global matrix 
system, where the subscripts indicate the matrix size required for N 
number of nodal points 



|^NxnJ j^Nxl| [^NxnJ |^Nxlj + J^NxNxN (NxN)xlj ^ 

2 

where ip is the time response of the flux at each nodal point and jp 

is the product ... N. Solved in this form, the major 

3 

limitation of the procedure was the large N size of the cubic array 

C, which is introduced as the result of the nonlinearity of Eq. (1). 

With this size requirement rapidly increasing for a finer mesh consisting 

3 

of more nodal points, the investigation of Reference 3 was limited to a 
small but practical N of 38. 



B. 



TREATMENT ON THE ELEMENT LEVEL 



The cubic form of C need not appear if it can be shown that the 
nonlinear portion of Eq. (8) can be broken apart into the product of 
two linear components such that one component can be evaluated on the 
element level and the element contributions then summed into global form, 



In general, for a nonlinearity of order m, a finite element 
problem of N nodes will generate a matirx size of N m+ 1. 
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resulting in a "regular” NxN matrix C* ready for multiplication by 
the remaining also summed to the global level. Proof that this is 



indeed the case, i.e., 



[^ C Nx(NxN)[[ V (NxN)xl} £ C NxnJ |%xl| 



is outlined as follows : 

Reference 3 established the element nonlinear term as 



- c yk ^ v 0 

A 

e 

where the integral is over the element area , 

C ijk = ff Wk (r e 4 e )drdz > i ' J ’ k>e = 1 - 2 ' 3 

A 

e 



( 9 ) 



(10) 



( 11 ) 



and £ being the local triangular coordinates, an innovation of Felippa 
[Ref. 6]. The coefficients C. ... form a 3x3x3 element array possessing 

ljk 

a regularity which suggests a kind of "cubic symmetry." The evaluation 
of this integral results in the following expression 

C lll ■ Te [ 24r l + 6 < r 2 + r 3 > ] 

C U2 ' Y<! [ 6r l + 4r 2 + 3r3 ] 

C 113 = Y t 6r l + 2r 2 + 4r 3 ] 

C 122 ‘ Y< “[ 4r l + 6l '2 + 2r 3 ] 

C 123 = Y [ 2 ^ r l + r 2 + r 3^1 
C 133 ' Y6 [ 4r l + 2r 2 + 6r 3l 
C 222 " YC [ 6r l + 24r 2 + 6r 3] 

C 333 ' Y6 [ 6r l + 6r 2 + 24r 3] 
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where 



and 



Y = ttA /180 
e 



r = r = r 
121 112 u 211 



C 321 C 231 C 213 etc as in Ref. 3. 



2 2 

Assuming that may be broken apart as = ip • <p 9 the LHS 

(left hand side) of Eq. (10) is rewritten as 

ff ^i<^j)^A )dA i,j,k = 1,2,3 



(13) 



or 



2tt 



Jf^ K ^1 C 2 C 3^ 



*1 




1 

i— l 


*2 


[ C 1 e 2 C 3 ] 


^2 


*3 




_*3_ 



r dr dz 



(14) 



Expanding and collecting terms yields 



2ir 



// 



’K 



h ri (r l C l + r 2 C l C 2 + r 3 ? l C 3^ 



+ ^2 (r l4^2 + W2 + 



+ ^3^ r l^l ? 3 + r 2 C 1^2 C 3 + r 3 C l C 3^' 



(15) 



+ ]^l' r l C l C 2 + r 2^1 C 2 + r 3 C l C 2 C 3^ + ^2^ r l ? l S 2 + r 2 C 2 + r 3^2 C 3^ 
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+ < f , 3^ r i ? i ? 2 ? 3 + r 2 ? 2 ? 3 + r 3 ? 2^3' > ( 

+ ^3 j ^ r l C 1^3 + r 2 C l C 2 C 3 + r 3 C l ? 3^ + lJ> 2^ r l ? l ? 2 ? 3 + r 2 C 2 C 3 + r 3^2^ 



2 2 3 

+ 4>3(r 1 C 1 C 3 + r 2 ^2 ? 3 + r 3^3^ 



dr dz 



k= 1, 2, 3 



Since the quantities ^ and <}>_. (i,j = 1,2,3) depend only on time, 
they are unaffected by the integration. Further, the expressions 
inside the parentheses are given a unique name d_^ to fl x their 
position in relation to a particular f° r some k = 1,2,3. Now 

the above expression may be written as 



2tt ff d. #1 dr dz . 

1 j JJ ijk 



(16) 



Using the integration formula from Felippa [Ref. 6] 



// 



£ m n £!ra!n! 

h h ? 3 dA ' Tm^+27! 



x 2A 



(17) 



where £,m,n are the exponents of £ in a specific d^_,^. The integration 
of (16) gives 



(ttA /180) ^ . (j) . f . .. 
e i j ljk 



(18) 



and f ^ x ^A^/ISO is identically equal to °f (12) 

The indicated summation is conveniently expressed by 



*3 f ijk = a i u ^ k = 1 » 2 » 3 



(19) 
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the desired 3x3 



e * 

which is immediately recognized as y ou ^ = 9 

element matrix. In this form it is combined into the global [CL, „] , 

NxN 

and the linearization has been accomplished. As will be shown later, 

it may be convenient to let <j> = when for small At, <f> = 

^t-At ~ ^t * This approximation is useful for predictor-corrector 

equation solvers in order to allow construction of [C., VT ] only once 

NxN 

for each successful time step. 
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III. SOLUTION PROCEDURE WITH THE LINEARIZED FORM 



For clarity, throughout the remainder of this thesis the following 

3 

naming convention will be used: N denotes the nonlinearized direct 

treatment by DVOGER as solved in Ref. 3. ’'Linearized approximate” 
indicates the linearized technique based on element level multiplication 
of C by 4> . . The "linearized exact" method, however, is the element 

level multiplication of C by the trial 4> t * The compact, or Nxp 
formulation may be used with any of the linearized methods in addition 
to the CRANKO treatment discussed later. 

A. GENERAL TREATMENT OF THE MATRIX EQUATION 

Solution of the matrix formulation of Eq. ( 8 ) involves clearing 

the LHS to provide a discrete relation for each 4^* This is accomplished 

either by iteration, which yields an approximation, or by matrix inversion, 

which has the advantage of giving the individual 4^ explicitly. In 

3 

both the nonlinearized N and the linearized NxN forms, Eq. ( 8 ) is 
reduced with multiplication by [A] thus 



[A] 1 [A] {ij>} = [A] 1 [B]{^} + [A] _ 1 [C^]{4.} 



( 20 ) 



in the linearized 4 > = 44 ; case or as 



[AfVm} = [A] 1 [B]{^} + [AfhCH* 2 } 



( 21 ) 



for the direct, nonlinear, formulation. This results in 



{*} = [AB]{i|i} + [A 



( 22 ) 
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where 



[C*] = [C*] or [C]W (23) 

as appropriate. The point here is that the [AB] and the [A are 
formed once and for all based on the time independent properties and 
geometry of the problem, whereas the [Cc}>] must be computed for each time 

t. Although the construction of the [C] is required only once in the 

3 2 

direct (N ) method, it must be additionally multiplied by {ty } at 

each time step. 

B. PRACTICAL CONSIDERATIONS WITH THE LINEARIZED TECHNIQUE 

Relating Eq. (22) to the demands of the subroutine DVOGER, which 
involves the evaluation of for every trial value of ip, it is 

apparent that the construction of [A] ^[C^]{^} for each of these trials 
could be a time consuming process, and hence a serious drawback of the 
linearized treatment. Noting, however, that for small At,ip ^ 
the approximation [Ci>] ^ [C^], where <j> = overcomes this 

disadvantage by allowing the construction of [A] ’*‘[Ccf>] only once for a 
given time t, regardless of hox^ many trial solutions are attempted by 
DVOGER. Equation (22) can now be further simplified since [A] ^[Ccj>] 
is readily combined with [AB] to reformulate the problem as 

• ^ ^ 

W = [c . ( 24 ) 

Thus, although the linearized form can be handled exactly, i.e., with 
<j> = tyy it may be more expeditous to establish the [C ] only once per 
valid time point, (that is, where the convergence criterion of the equation 
solution has been satisfied). In this case, the linearized treatment is 
then an approximate technique in that <$> = ^ 
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c. 



DIFFICULTIES WITH THE APPROXIMATION 



The simplification expressed as Eq. (24) works extremely well when 
used with DVOGER as long as large changes in ip are incurred by small 

At (the so-called "stiff” system). When the solution approaches steady 

# 

state, however, ip approaches zero faster than (^ - ^t-At^^ t does, 
and this causes a second kind of perturbation wherein the system becomes 

increasingly sensitive to small changes in ij; - Tp . Apparently 

• • 

DVOGER requires very accurate information. The ip(t) provided from 

ijj(t) x does not match the DVOGER prediction of which is 

• 2 

based on the assumption that was generated by This 

discrepancy in the neighborhood of ip 0 causes DVOGER to be unable to 
recognize the steady state solution. Prepared especially for stiff 
systems, the equation-solver expects, instead, a large scale change in 
ip. The time step is therefore reduced accordingly, which means that very 
slow progress is made in the steady-state region which is otherwise 
handled quite rapidly under the exact methods. The implication here is 
that, in those cases which reach terminal flux values rather early in 
problem time, the exact formulation of the linearized treatment is to be 



preferred, even though the [C ] must be computed for each attempted 



. This aspect has been analyzed for the test problems considered, 
with the results given in Table II. 
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IV. THE Nxp COMPACTING SCHEME 



A. ESTABLISHED COMPACTING METHOD 

The problem under consideration here, being of the form A_^^ = 

B (i|>i ... i|^) , is formally handled by multiplying through by [A] \ 
discussed earlier. For large matrices, the inversion process, and indeed 
all matrix operations, become extremely costly in terms of machine time 
and storage. The fact that the matrices may be banded or symmetric does 
offer significant economy when the matrix is not inverted. For example, 
a symmetric NxN matrix may be stored as NxN/2. For matrices resulting 
from finite element formulation, the numbering schedule determines a band 
width to be used with conventional compacting schemes. In the finite 
element case, the difference is compared for each node in the mesh. The 
largest of these differences plus one equals the minimum band width 
allowed. For example, node 7 in Figure 1 has a maximum difference of 
13 - 2 = 11, whereas node 22 has the largest difference (37 - 10 = 27) 
for the system. Thus the resultant minimum band width would be 28. As 
the mesh increases in nodal points, the minimum band width must also 
increase in response, no matter how adroitly the numbers are assigned. 
Regardless of how sparce a matrix may be, these schemes all ensure that 
the "compactness" grows commensurately . In addition, the calculation of 
an inverse becomes ever more unpleasant, especially since this operation 
on a sparce matrix generally results in a dense, nonsymmetric inverse. 
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B. BASIS FOR THE COMPACTING SCHEME 

In solving a matrix equation A^i^ = ••• 5 it is not 

necessary to form the inverse, but rather, multiply out the left hand 
side and rearrange to obtain a set of equations: 



ip. = (B.-A. .*.)/A. . 

i l lj j li 



(25) 



where i,j,k =1, . .. N. This procedure is especially suitable for 
sparce systems since the number of non-zero A^ entities will be small. 
Further, if these non-zero A^ can be located and tagged with an 
identifier, then they may be stored together in a dense array, thereby 
replacing the standard NxN size matrix by an equivalent array of size 
Nxp, where p is the maximum number of non-zero entries on any one row. 

This type of compacting is especially well suited to matrices 
resulting from finite element formulation, with p determined by the 
choice of the finite element and the discretized model. In the case 
under study, a linear triangular element is used such that, when assembled 
on the global level, it forms, for the purpose of illustration, a network 
of interlocking hexagonal polygons, each having an apex located directly 
over its parent central node. Such a pattern is sketched as Figure 2. 

Note that the sides are indeed the contributing "neighbor” nodes. Boundary 
nodes will, of course, be missing any "would be" contributors sought in 
the boundary exterior. In a regularly drawn system, as in Figure 2, an 
interior node will have six neighbors plus itself for a total of seven 
contributors. That is to say, that no matter how large the system of 
mesh points, a relation describing the effect of the system on any given 
point would contain a maximum of seven non-zero values. 
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Geometric considerations or a desire to examine some area of the 



mesh more closely may result in systems of irregular discretization 
wherein the maximum number of entries may be greater than seven. This 
analysis may be extended to consider other finite element shapes in a 
similar manner. 

C. CONSTRUCTION OF THE Nxp ARRAY 

To utilize this attribute of finite elements to produce the dense 
Nxp system, it is necessary only to construct a table of nodal points 
listing their contributors. This table, or connectivity array, is then 
used as an index to locate the values of the various quantities associa- 
ted with a particular node. All the standard matrix operations can be 
performed on these compacted arrays, but the construction of a matrix 
inverse has no usefulness since a NxN system is then reconstituted from 
an Nxp array. 

As an example, construction of the connectivity vector associated 
with node 18 in Figure 1 is 

[18 17 27 26 19 14 13 12] 

and the vector associated with node 6 is 

[6 11 7 21 0 0 0 0]. 

The contributors may be entered in any order except that, for computa- 
tional ease, the central, or "parent" node is the first element. These 
vectors are collected into the node neighbor connectivity array of size 
Nxp, or, for Figure 1, 38x8. 
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To multiply a Nxp array of a_^ elements by a vector of 
elements, a search of the connectivity array is performed on the "i"th 
row to locate a match for the b^ . When the match is found, the values 
represented by a_^ and b^ are multiplied together resulting in a new 
CL vector. If the system is correctly formed and compacted, there 
always will be a match. 

The compacting scheme is bound by the geometry of the problem, and 
as such, has meaning only as it pertains to that problem. For example, 
operations with two equi-sized Nxp arrays of different connectivity 
relations cannot be performed. This is in contrast to regular matrix 
operations where the restriction is only to size and not to origin. A 
short computer program to demonstrate the operation of this Nxp scheme 
is given, with results, in Appendix C. 
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V. THE CRANK-NICOLSON FORMULATION 



Recognizing the very sparce nature of the system of differential 
equations formulated by this problem and, if the matrix inversion process 
is to be avoided, the necessity of iteration, a straightforward attack 
based on the definition of the derivative (finite difference) appears as 
a feasible alternative to the DVOGER technique. One such approach was 
presented by J. Crank and P. Nicolson in 1947 [Ref. 7] and has been 
discussed in many works throughout the ensuing years. This method has 
been shown to be unconditionally stable. An adaptation of this method 
is given here as follows: 

In the general sense, a matrix formulation of a set of linear 
differential equations, can be represented as 



where ^ is a function of time. For the case of an initial disturbance 
as the forcing function, F(t) = 0 for t > 0. Focusing attention for 
the moment on a particular equation and expressing [16] in terms of the 
definition of the derivative 




• • • , 



N 



(26) 





(27) 



which is 



[A- 1/2* At* = [A + (1/2* At* C) 



(28) 



[(2/AfA)-C]tJ» t = [(2/At*A)+C]^ t _ At 



(29) 
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letting D = (2/At*A)-C, E = (2/At*A)+C, <j> = ^t-At t * ien > returning 

to the system form the result is 



D. = E . . <}> . = E.* 

ij ij Y j i 



i> j = 1, 



N . 



(30) 



This constitutes a set of linear simultaneous equations which may be 
solved by any convenient method. Chosen here is Gauss-Seidel iteration 
[Ref. 8] which is well suited for the problem at hand since, the sum of 

and E_^_. <f>^. contain very few non-zero terms, thereby incurring 

relatively small roundoff errors regardless of system size. 

This procedure, which also uses the Nxp compacting, has been 
written for this thesis as the computer code CRANKO. After receiving 
the matrix information and control parameters from the calling program, 
CRANKO first builds the C matrix of Eq. (26) for the initial cf> value. 
Then it selects a trial time interval and, forming the equation set (30) 
based on this At, attempts to find a solution satisfying a specified 
convergence criterion based on a relative error. If no satisfactory 
solution set is found after a number of iterations, a smaller At is 
selected; and, after recomputation of Eq. (30), another attempt at 
convergence is made. The successful solutions then replace the previous 
ip thereby forming the starting point for another cycle, beginning with 
a new computation of the C matrix. The choice of At is the control- 
ling factor. The scheme employed counts the number of iterations required 
for a solution. When an increasing number of iterations indicate greater 

i 

difficulty, the interval At is decreased; whereas, if fewer iterations 
are required, the time step may be relaxed and made larger. Experience 
with the CRANKO routine establishes this method as an extremely fast 
technique for solving this type of problem. 
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For those readers unfamiliar with the Gauss-Seidel iteration scheme, 
a simple example using a 3x3 system is given: 



Vi = E j * ■ Vat 



d ll d l2 d 13 




1 

i— 1 

1 




I— 1 

CD 

1 


d 21 d 22 d 23 




CM 

-a- 


tl 


e 2 


— d 31 d 32 d 33- 




— i 

Ol 

-2- 

I 




- e 3- 



first iteration: 



d ll*l + d 12*2 + d 13*3 = e l 



= [e, - 



= 



( d 12 < f>2 + d i3*3^/ d n 



d 21*l + d 22^2 + d 23^3 



= e 



2 

[e. 



( d 2i'J ; i + d 23^3^^ d 22 



= 4 >< 



d 31*l + d 32^2 + d 33^3 



= e. 



*3 ‘ [e 3 



(d 31*l + d 32*2 )1/d 33 



for successive iterations, 
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d ll*l + d 12*2 + d 13*3 = e l 



^1 ^ e l ^ d 12^2 + d 13^3' ) ^ d ll 
** 

*1 ‘ *1 



d 21*l + d 22^2 + d 13*3 ' e ; 



*2 ° |e 2 ‘ (d 21+l + d 13 ,| '3 ),/d 22 



4 '**_ 

2 " 2 



ft £ 

d 31*l + d 32*2 + d 33*3 " e 3 



*3 ’ [e 3 ‘ <d 31*l + d 32*2 >1/d 33 



etc. 
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VI. SUMMARY OF AVAILABLE METHODS 



For clarity, the innovative computational tools introduced in this 

thesis are summarized. First is the linearizing of the governing equation 

Eq, (8) by the multiplication on the element level of by ^ . 

From this development, two options appear feasible: either the exact 

* 

computation of [C ] for each trial ip, or the approximation of 

* 



ip zz to allow the construction of [C ] only once during the 

search for a new 1 J 1 . This approach uses, as does the "direct," or 
nonlinear treatment involving [C], the inverting of [A] to clear the 
LHS as shown in Eq, (22), 

The second innovation is the reduction of the NxN matrix to a 
smaller Nxp compact form based, not upon the mathematics of the problem, 
but upon the geometry of the finite element model. The same solution 
methods may be used, with the only difference being that the inverting 
of [A] is replaced by iteration. This iteration can only be as accurate 
as the imposed error criterion. For large systems, however, where the 
[A] ^ produces a full NxN matrix on the RHS, this iteration is 
especially beneficial. 

The last contribution was the development of a new equation-solver 
based on a different theory than the DVOGER routine. The new code deals 
with the differential equations directly and employs the Nxp compacting. 
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VII. TEST PROBLEMS AND RESULTS 



A. THE PHYSICAL MODEL 

A cylindrical nuclear reactor consisting of a core region and 
blanket having overall a height of 220 cm and a diameter of 180 cm was 
modeled for the test problems. As shown in Figure 1, a radial slice 
from this reactor was discretized by a finite element of 54 elements and 
38 nodal points. After a detailed analysis of this mesh was completed, 
preliminary investigation was begun using a model of 220 elements 
(Figure 3). 

B. COMPUTER PROCESSING CONSIDERATIONS 

All programs were written in the FORTRAN IV language and processed 
on the IBM 360/67 computer using the FORTRAN f H f compiler. The results 
presented in terms of CPU storage requirements and processing time ought 
to be indicative, in a relative way, of those obtained from other machines. 
Single precision (seven significant digits) was used throughout the 
investigation, with the results later verified using double precision 
(fifteen significant digits). The observed relative discrepancy was less 
than .01%. It is recongized, however, that for larger systems a more 
significant difference may result. 

C. PROBLEM ANALYSIS 

The three techniques discussed have been incorporated into computer 

codes and have been tested on those problems described in Reference 3, 

namely, the dynamic reactor response to 1) a disturbance at the center, 

2) a uniform disturbance throughout the core, and 3) a disturbance occurring 

3 

at a skew point in the core (R - 40 cm, Z ~ 0 cm). The direct, or N 
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matrix (an Nx(N ) array) equation was solved in Reference 3, and those 

results are used in this thesis as a standard for comparison. The 

computer codes developed here have been constructed within the framework 

3 

and parameters of the original N program, thereby allowing any 
discrepancy to be attributed to the technique employed. For each tech- 
nique and disturbance, the flux at three sample points was recorded 
throughout the time history of the solutions and written onto a computer 
storage device. A separate program was written to process this infor- 
mation and produce the correlations in tabular and graphic form. The 
graphic results are presented in Figures 4-21 for the center point 
(R = 0 cm, Z = 0 cm), a core point (R = 40 cm, Z = 40 cm), and a reflector 
point (R = 75 cm, Z = 80 cm). Table II lists the specific programs 
employed arranged by disturbance type and gives parameters such as 
storage requirements and computer processing time. 

1. Solutions with the CRANKO Subroutine 

All techniques appear to give acceptable accuracy, however 
there is a startling difference in processing time. The CRANKO subroutine 
was able to accomplish in seconds the work which required minutes for the 
DVOGER subroutine (See Table II). The steady state solution given by 
CRANKO was in every case higher than that obtained by DVOGER. The 
implication here is not clear since there is no analytical solution 
available to represent the exact answer. In general, the transient 
solutions agreed well, except for the oscillation of the CRANKO solution 
in the early stages, most clearly depicted in the extreme case of Figures 
8 and 9. This fluctuation is most likely the result of inadequate control 
of the CRANKO time step selection where apparently an interval has been 
chosen which is too large to maintain good accuracy. It is noteworthy, 
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however, that the technique is so stable that even after large deviation 
it is able to recognize the error, correct itself, and return to the 
solution curve given by the DVOGER routine. 

2. Solutions with the DVOGER Subroutine 

The linearized technique using premultiplication by the previous 
time step flux conformed to theoretical expectations. As the 

steady state solution was approached, this method proceeded with in- 
creasing difficulty. Hence much more processing time to advance through 

the steady state was required. 

3 

The direct N method, on the other hand, has no difficulty 

in recognizing the steady state and marches forward with large time 

increments until problem termination. A very successful effort was made 

to use the exact treatment with the linearized technique. As expected, 

more processing time was consumed in getting to the steady state due to 

•k 

the requirement to recompute the [C for each trial ifj, but once the 

steady state was approached, the procedure progressed rapidly as in the 
direct treatment. Comparing the direct with the linearized exact method, 
theory developed in this thesis predicts identical results with twice 
the computing time required by the linearization. As shown in Table II, 
this is indeed the case. The largest relative error observed between the 
two methods was .06%, which is attributed to roundoff steming from the 
doubling of the computational requirements. When the greatest accuracy 
is required, the exact treatment with linearization may therefore prove 
the most useful, especially if the problem rapidly reaches steady state. 
This discussion of the linearized technique is applicable both to the NxN 
method or to the Nxp compact method. However since the Nxp method 
requires iteration, more processing time and less accuracy might be 
expected. 
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D. EVALUATION OF ERROR 

The error analysis for these procedures must remain subjective as 
there is no "correct 11 result. Several sources of error are present and 
deserve comment. Tests were made with both EPS = .1 and EPS = .01, the 
error criterion for the DVOGER routine, and little difference was noted 

for the large increase in processing time. [EPS is defined by DVOGER as 

r ERROR (i) n r , . , . . 

I — m a J. The convergence criterion [relative deviation of any 

i=l MAX 

^ between successive iterations] for the iterative solutions in CRANKO 

was varied by tenths from .01 to .001 with the result of producing higher 

steady state values of only a few percent as the criterion was made more 

stringent. Processing time for this variation increased only a few 

seconds. Double percision trials for the linearized DVOGER method and 

the compact CRANKO were made with neglibible difference. 

The grid of only 54 elements theoretically introduces the largest 

source of error from the "true" solution. With the techniques developed 

in this thesis, a finer grid size of 220 elements has been implemented 

/ 

into the CRANKO and the linearized codes. Because the linearized code 
(DVOGER) requires about 566K core size, extensive testing with this 
technique is not anticipated, but the CRANKO routine which requires only 
188K is currently being investigated. As the results of this investiga- 
tion become clear, a better understanding of the error induced by each of 
the various approaches to this problem would be obtained. 
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VIII. CONCLUSION 



It appears that nonlinear differential equations resulting from a 
finite element formulation (or for that matter finite differencing) may 
be succesfully reduced using the linearizing technique developed in 
this thesis. The compacting scheme for matrix equations effects gross 
savings for storage of large systems (Table I) and becomes even more 
attractive when the inverse of a matrix is either not desired or required. 
The Crank-Nicolson formulation coupled with Gauss-Seidel iteration as 
expressed in the computer code CRANKO has demonstrated itself to be a 
viable, extremely fast technique for the solution of differential 
equations. Further improvements to the CRANKO routine, specifically in 
regard to better control of the time step selection, may improve upon the 
attractiveness of this approach. 

Development is warranted in several areas. The computation of the 
Jacobian using iteration must be formulated in order to use the DVOGER 
subroutine with the compact system. DVOGER itself might usefully be 
reviewed to reduce its size requirements also. Investigation of the 
problems using the fine mesh size of 220 elements is continuing. This 
should yield more exact results and offer a better standard for comparison. 
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APPENDIX A 



TABLES 



TABLE I 



IBM 360/67 CORE REQUIREMENTS 
in K Bytes 



SYSTEM 


NO. NODES 


OVERHEAD 


MATRIX STORAGE 


TOTAL 


NxN 2 


38 


100 


225 


325 


NxN 


38 


100 


35 


135 


Nxp 


38 


100 


12 


112 


NxN 


132 


150 


418 


568 


Nxp 


132 


150 


38 


188 


Nxp 


400 


290 


112 


402 t 



-(•calculated estimate 
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